Overview

Data

Meta-ecosystems

All

Equally sized

Biomass

p_boxplots = ds_metaecosystems %>%
  mutate(metaecosystem_type = fct_reorder(metaecosystem_type, total_metaecosystem_bioarea)) %>%
  filter(disturbance == disturbance_input,
         metaecosystem_type %in% c("S_S", "M_M", "L_L")) %>%
  ggplot (
    aes(
      x = day,
      y = total_metaecosystem_bioarea,
      group = interaction(day, metaecosystem_type),
      fill = metaecosystem_type
    )
  ) +
  geom_boxplot() +
  labs(x = "Day",
       y = "Meta-ecosystem Total Bioarea (µm²)",
       color = '',
       fill = '') +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small-Small",
               "Medium-Medium",
               "Large-Large")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )  +
  scale_x_continuous(breaks = unique(ds_metaecosystems$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(caption = "Vertical grey line: resource flow")

p_boxplots

Beta Diversity

p_boxplots = ds_metaecosystems %>%
  mutate(metaecosystem_type = fct_reorder(metaecosystem_type, total_metaecosystem_bioarea)) %>%
  filter(disturbance == disturbance_input,
         metaecosystem_type %in% c("S_S", "M_M", "L_L")) %>%
  ggplot (aes(
    x = day,
    y = bray_curtis,
    group = interaction(day, metaecosystem_type),
    fill = metaecosystem_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) + 
  labs(
    x = "Day",
    y = "Beta diversity (Bray-Curtis Index)",
    color = '',
    fill = '',
    caption = "Vertical grey line: resource flow"
  ) +
   scale_fill_manual(
     values = c(colour_small, colour_medium, colour_large),
     labels = c("Small-Small",
                "Medium-Medium",
                "Large-Large")
   ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

p_boxplots

Gamma Diversity

p_boxplots = ds_metaecosystems %>%
  mutate(metaecosystem_type = fct_reorder(metaecosystem_type, total_metaecosystem_bioarea)) %>%
  filter(disturbance == disturbance_input,
         metaecosystem_type %in% c("S_S", "M_M", "L_L")) %>%
  ggplot (aes(
    x = day,
    y = gamma_shannon,
    group = interaction(day, metaecosystem_type),
    fill = metaecosystem_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) + 
  labs(
    x = "Day",
    y = "Gamma diversity (Shannon Index)",
    color = '',
    fill = '',
    caption = "Vertical grey line: resource flow"
  ) +
   scale_fill_manual(
     values = c(colour_small, colour_medium, colour_large),
     labels = c("Small-Small",
                "Medium-Medium",
                "Large-Large")
   ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

p_boxplots

MM vs SL

Biomass

ds_metaecosystems %>%
      filter (disturbance == disturbance_input,
              metaecosystem_type %in% c("S_L", "M_M")) %>%
      ggplot (
        aes(
          x = day,
          y = total_metaecosystem_bioarea,
          group = system_nr,
          fill = system_nr,
          color = system_nr,
          linetype = metaecosystem_type
        )
      ) +
      geom_line () +
      labs(
        x = "Day",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        fill = "System nr",
        color = "System nr",
        linetype = ""
      ) +
      scale_linetype_discrete(labels = c("Medium-Medium",
                                         "Small-Large")) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_x_continuous(breaks = unique(ds_metaecosystems$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      labs(caption = "Vertical grey line: resource flow")

p_boxplots = ds_metaecosystems %>%
      filter(disturbance == disturbance_input,
             metaecosystem_type %in% c("S_L", "M_M")) %>%
      ggplot (
        aes(
          x = day,
          y = total_metaecosystem_bioarea,
          group = interaction(day, metaecosystem_type),
          fill = metaecosystem_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      )  +
      scale_x_continuous(breaks = unique(ds_metaecosystems$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      labs(caption = "Vertical grey line: resource flow")

p_boxplots

filtered_data = ds_metaecosystems %>%
    filter(metaecosystem_type %in% c("M_M", "S_L"),
           disturbance == disturbance_input,
           time_point >= first_time_point_model,
           time_point <= last_time_point_model)
full_model = lmer(
  total_metaecosystem_bioarea ~
    day +
    metaecosystem_type +
    metaecosystem_type : day +
    (day | system_nr),
  data = filtered_data,
  REML = FALSE,
  control = lmerControl(optimizer = "Nelder_Mead")
)

plot(full_model)

qqnorm(resid(full_model))

null_model = lmer(
  total_metaecosystem_bioarea ~
    day +
    (day | system_nr),
  data = filtered_data,
  REML = FALSE,
  control = lmerControl(optimizer = "Nelder_Mead")
)

AIC = AIC(full_model, null_model) #When I add baseline_bioarea : day, the AIC of the model stays the same when I take off the ecosystem type.

AIC
##            df      AIC
## full_model  8 1440.690
## null_model  6 1439.892
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Data: filtered_data
## Models:
## null_model: total_metaecosystem_bioarea ~ day + (day | system_nr)
## full_model: total_metaecosystem_bioarea ~ day + metaecosystem_type + metaecosystem_type:day + (day | system_nr)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## null_model    6 1439.9 1452.5 -713.95   1427.9                     
## full_model    8 1440.7 1457.4 -712.35   1424.7 3.2021  2     0.2017
p_value = anova$`Pr(>Chisq)`[2]
p_value = round(p_value, digits = 3)
p_boxplots + 
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value))

Time point 1

Let’s start by looking at whether meta-ecosystem type and disturbance had an effect at time point = 1. At this time point, no disturbance event had yet occurred. Therefore, we would not expect an effect of disturbance. In regards to meta-ecosystem type, there might be an effect if it comes from just the sizes of the two ecosystems.

chosen_time_point = 1

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_1 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_1

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 243.5524
## null_model  2 242.3665
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1      8 1.2142e+10                             
## 2      9 1.3172e+10 -1 -1029927000 0.6786 0.4339
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_1 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Time point 2
chosen_time_point = 2

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_2 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_2

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 245.2379
## null_model  2 243.3101
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df  Sum of Sq     F Pr(>F)
## 1      8 1.4371e+10                           
## 2      9 1.4475e+10 -1 -104202831 0.058 0.8157
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_2 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Time point 3
chosen_time_point = 3

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_3 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_3

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 240.6479
## null_model  2 247.4085
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df   Sum of Sq      F  Pr(>F)  
## 1      8 9.0814e+09                                
## 2      9 2.1808e+10 -1 -1.2727e+10 11.211 0.01011 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_3 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Time point 4
chosen_time_point = 4

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_4 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_4

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 248.7535
## null_model  2 247.2519
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1      8 2.0426e+10                             
## 2      9 2.1469e+10 -1 -1043835025 0.4088 0.5404
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_4 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Time point 5
chosen_time_point = 5

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_5 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_5

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 240.9001
## null_model  2 239.1289
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1      8 9313404235                            
## 2      9 9528897119 -1 -215492884 0.1851 0.6784
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_5 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Time point 6
chosen_time_point = 6

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_6 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_6

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 226.1687
## null_model  2 224.6854
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1      8 2134661375                            
## 2      9 2247857851 -1 -113196476 0.4242 0.5331
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_6 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Time point 7
chosen_time_point = 7

filtered_data = ds_metaecosystems %>%
                            filter(disturbance == disturbance_input,
                                   metaecosystem_type %in% c("M_M", "S_L"),
                                   time_point == chosen_time_point)

day_of_time_point = unique(filtered_data$day)
time_point_7 = filtered_data %>%
      ggplot (aes(x = metaecosystem_type,
                  y = total_metaecosystem_bioarea)) +
      geom_boxplot() +
      labs(
        title = paste0(
          "Time point = ",
          chosen_time_point
        ),
        x = "",
        y = "Meta-ecosystem Total Bioarea (µm²)",
        color = '',
        fill = ''
      ) +
      scale_fill_manual(
        values = c(colour_MM, colour_SL),
        labels = c("Medium-Medium",
                   "Small-Large")
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

time_point_7

full_model = lm(total_metaecosystem_bioarea ~
            metaecosystem_type,
          data = filtered_data)

par(mfrow=c(2,3))
plot(full_model, which = 1:5)

null_model = lm(total_metaecosystem_bioarea ~
                  1,
          data = filtered_data)


AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  3 219.4978
## null_model  2 218.0579
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: total_metaecosystem_bioarea ~ metaecosystem_type
## Model 2: total_metaecosystem_bioarea ~ 1
##   Res.Df        RSS Df Sum of Sq      F Pr(>F)
## 1      8 1095508451                           
## 2      9 1158623185 -1 -63114734 0.4609 0.5164
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

time_point_7 +
  labs(
    title = paste0(
      "Day = ",
      day_of_time_point,
      ", ΔAIC = ",
      deltaAIC ,
      ", p = ",
      p_value,
      ", Adjusted R2 = ",
      R2_P
    )
  )

Beta Diversity

ds_metaecosystems %>%
  filter (disturbance == disturbance_input,
          metaecosystem_type %in% c("S_L", "M_M")) %>%
  
  ggplot (
    aes(
      x = day,
      y = bray_curtis,
      group = system_nr,
      fill = system_nr,
      color = system_nr,
      linetype = metaecosystem_type
    )
  ) +
  geom_line () +
  labs(
    title = paste("Disturbance =", disturbance_input),
    x = "Day",
    y = "Beta Diversity (Bray-Curtis Index)",
    fill = "System nr",
    color = "System nr",
    linetype = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_linetype_discrete(labels = c("Medium-Medium",
                                     "Small-Large")) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_x_continuous(breaks = unique(ds_metaecosystems$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  )

p_boxplots = ds_metaecosystems %>%
  filter(disturbance == disturbance_input,
         metaecosystem_type %in% c("S_L", "M_M")) %>%
  ggplot (aes(
    x = day,
    y = bray_curtis,
    group = interaction(day, metaecosystem_type),
    fill = metaecosystem_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) + 
  labs(
    x = "Day",
    y = "Beta diversity (Bray-Curtis Index)",
    color = '',
    fill = '',
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_MM, colour_SL),
    labels = c("Medium-Medium",
               "Small-Large")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

p_boxplots

filtered_data = ds_metaecosystems %>%
    filter(time_point >= 2,
           disturbance == disturbance_input,
           metaecosystem_type %in% c("M_M", "S_L"))
full_model = lmer(
  bray_curtis ~
    day +
    metaecosystem_type +
    metaecosystem_type : day +
    (day | system_nr),
  data = filtered_data,
  REML = FALSE,
  control = lmerControl(optimizer = "Nelder_Mead")
)

summary(full_model)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bray_curtis ~ day + metaecosystem_type + metaecosystem_type:day +  
##     (day | system_nr)
##    Data: filtered_data
## Control: lmerControl(optimizer = "Nelder_Mead")
## 
##      AIC      BIC   logLik deviance df.resid 
##    -95.4    -78.7     55.7   -111.4       52 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.31770 -0.80339 -0.00421  0.63375  2.23000 
## 
## Random effects:
##  Groups    Name        Variance  Std.Dev.  Corr 
##  system_nr (Intercept) 1.658e-12 1.288e-06      
##            day         6.670e-15 8.167e-08 -1.00
##  Residual              9.138e-03 9.559e-02      
## Number of obs: 60, groups:  system_nr, 10
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                0.053272   0.049188   1.083
## day                        0.008414   0.002555   3.293
## metaecosystem_typeS_L     -0.029385   0.069563  -0.422
## day:metaecosystem_typeS_L  0.006708   0.003613   1.857
## 
## Correlation of Fixed Effects:
##             (Intr) day    mt_S_L
## day         -0.935              
## mtcsyst_S_L -0.707  0.661       
## dy:mtcs_S_L  0.661 -0.707 -0.935
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
plot(full_model)

qqnorm(resid(full_model))

null_model = lmer(
  bray_curtis ~
    day +
    (day | system_nr),
  data = filtered_data,
  REML = FALSE,
  control = lmerControl(optimizer = "Nelder_Mead")
)

AIC = AIC(full_model, null_model) 
AIC
##            df       AIC
## full_model  8 -95.44362
## null_model  6 -85.84173
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Data: filtered_data
## Models:
## null_model: bray_curtis ~ day + (day | system_nr)
## full_model: bray_curtis ~ day + metaecosystem_type + metaecosystem_type:day + (day | system_nr)
##            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## null_model    6 -85.842 -73.276 48.921  -97.842                        
## full_model    8 -95.444 -78.689 55.722 -111.444 13.602  2   0.001113 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_value = anova$`Pr(>Chisq)`[2]
p_value = round(p_value, digits = 3)
MM_SL_beta_diversity = p_boxplots + 
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value))

MM_SL_beta_diversity

Gamma Diversity

ds_metaecosystems %>%
  filter (disturbance == disturbance_input,
          metaecosystem_type %in% c("S_L", "M_M")) %>%
  
  ggplot (
    aes(
      x = day,
      y = gamma_shannon,
      group = system_nr,
      fill = system_nr,
      color = system_nr,
      linetype = metaecosystem_type
    )
  ) +
  geom_line () +
  labs(
    title = paste("Disturbance =", disturbance_input),
    x = "Day",
    y = "Beta Diversity (Bray-Curtis Index)",
    fill = "System nr",
    color = "System nr",
    linetype = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_linetype_discrete(labels = c("Medium-Medium",
                                     "Small-Large")) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_x_continuous(breaks = unique(ds_metaecosystems$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  )

p_boxplots = ds_metaecosystems %>%
  filter(disturbance == disturbance_input,
         metaecosystem_type %in% c("S_L", "M_M")) %>%
  ggplot (aes(
    x = day,
    y = gamma_shannon,
    group = interaction(day, metaecosystem_type),
    fill = metaecosystem_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Gamma diversity (Shannon Index)",
    color = '',
    fill = '',
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_MM, colour_SL),
    labels = c("Medium-Medium",
               "Small-Large")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )
  
p_boxplots

filtered_data = ds_metaecosystems %>%
    filter(time_point >= 2,
           disturbance == disturbance_input,
           metaecosystem_type %in% c("M_M", "S_L"))
full_model = lmer(
  gamma_shannon ~
    day +
    metaecosystem_type +
    metaecosystem_type : day +
    (day | system_nr),
  data = filtered_data,
  REML = FALSE,
  control = lmerControl(optimizer = "Nelder_Mead")
)

summary(full_model)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: gamma_shannon ~ day + metaecosystem_type + metaecosystem_type:day +  
##     (day | system_nr)
##    Data: filtered_data
## Control: lmerControl(optimizer = "Nelder_Mead")
## 
##      AIC      BIC   logLik deviance df.resid 
##      1.8     18.5      7.1    -14.2       52 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2502 -0.5035  0.0383  0.7035  1.7116 
## 
## Random effects:
##  Groups    Name        Variance  Std.Dev. Corr 
##  system_nr (Intercept) 1.646e-02 0.128291      
##            day         6.973e-05 0.008351 -1.00
##  Residual              4.308e-02 0.207551      
## Number of obs: 60, groups:  system_nr, 10
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                1.568064   0.121231  12.935
## day                       -0.004521   0.006687  -0.676
## metaecosystem_typeS_L      0.116940   0.171447   0.682
## day:metaecosystem_typeS_L -0.001375   0.009457  -0.145
## 
## Correlation of Fixed Effects:
##             (Intr) day    mt_S_L
## day         -0.948              
## mtcsyst_S_L -0.707  0.670       
## dy:mtcs_S_L  0.670 -0.707 -0.948
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
plot(full_model)

qqnorm(resid(full_model))

null_model = lmer(
  gamma_shannon ~
    day +
    (day | system_nr),
  data = filtered_data,
  REML = FALSE,
  control = lmerControl(optimizer = "Nelder_Mead")
)

AIC = AIC(full_model, null_model)
AIC
##            df       AIC
## full_model  8 1.7795503
## null_model  6 0.5805223
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Data: filtered_data
## Models:
## null_model: gamma_shannon ~ day + (day | system_nr)
## full_model: gamma_shannon ~ day + metaecosystem_type + metaecosystem_type:day + (day | system_nr)
##            npar     AIC    BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model    6 0.58052 13.147 5.7097  -11.419                    
## full_model    8 1.77955 18.534 7.1102  -14.220 2.801  2     0.2465
p_value = anova$`Pr(>Chisq)`[2]
p_value = round(p_value, digits = 3)
MM_SL_gamma_diversity = p_boxplots + 
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value))

MM_SL_gamma_diversity

SL vs SL*

Patches

All

Isolated

Biomass

ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = bioarea_per_volume,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Patch Bioarea (µm²/μl)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Alpha Diversity

ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = shannon,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Alpha Diversity (Shannon Index)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Species Densities

Ble
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Ble,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Ble Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Cep
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Cep,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Cep Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Col
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Col,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Col Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Eug
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Eug,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Eug Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Eup
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Eup,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Eup Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Lox
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Lox,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Lox Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Pau
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Pau,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Pau Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Pca
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Pca,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Pca Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Spi
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Spi,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Spi Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Spi_te
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Spi_te,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Spi_te Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Tet
ds_patches %>%
  filter(disturbance == disturbance_input,
       metaecosystem == "no") %>%
  mutate(eco_metaeco_type = fct_reorder(eco_metaeco_type, patch_size_volume)) %>%
  ggplot(aes(
    x = day,
    y = Tet,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Tet Density (indiv/μl ?)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_fill_manual(
    values = c(colour_small, colour_medium, colour_large),
    labels = c("Small isolated",
               "Medium isolated",
               "Large isolated")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Small

Biomass

ds_patches %>%
      filter(disturbance == disturbance_input,
             patch_size == "S") %>%
      ggplot(
        aes(
          x = day,
          y = bioarea_per_volume,
          group = culture_ID,
          fill = culture_ID,
          color = culture_ID,
          linetype = eco_metaeco_type
        )
      ) +
      geom_line(stat = "summary", fun = "mean") +
      labs(
        x = "Day",
        y = "Patch Bioarea (µm²/μl)",
        linetype = ""
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_linetype_discrete(
        labels = c(
          "Small isolated",
          "Small connected to small",
          "Small connected to large"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      labs(caption = "Vertical grey line: resource flow")

ds_patches %>%
      filter(patch_size == "S",
             disturbance == disturbance_input) %>%
      ggplot(
        aes(
          x = day,
          y = bioarea_per_volume,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Patch Bioarea (µm²/μl)",
        fill = ""
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      ) +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      labs(caption = "Vertical grey line: resource flow")

p_effect_size = ds_patches_effect_size %>%
  filter(
    !time_point == 0,
    disturbance == disturbance_input,
    eco_metaeco_type %in% c("S (S_S)", "S (S_L)")
  ) %>%
  ggplot(aes(x = day,
             y = bioarea_per_volume_d,
             color = eco_metaeco_type)) +
  geom_point(position = position_dodge(0.5)) +
  geom_line(position = position_dodge(0.5)) +
  geom_errorbar(
    aes(ymin = bioarea_per_volume_d_lower,
        ymax = bioarea_per_volume_d_upper),
    width = .2,
    position = position_dodge(0.5)
  ) +
  labs(
    x = "Day",
    y = "Bioarea Density Effect Size (Hedge's d)",
    color = ""
  ) +
  scale_color_discrete(labels = c("Small connected to large",
                                  "Small connnected to small")) +
  scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted",
    color = "black",
    size = 0.7
  )

p_effect_size

filtered_data = ds_patches_effect_size %>%
                         filter(time_point >= first_time_point_model,
                                time_point <= last_time_point_model,
                                disturbance == disturbance_input,
                                eco_metaeco_type %in% c("S (S_S)","S (S_L)"))
full_model = lm(bioarea_per_volume_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(bioarea_per_volume_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 32.88298
## null_model  2 48.33802
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: bioarea_per_volume_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: bioarea_per_volume_d ~ 1
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
## 1      8  4.7301                                
## 2     11 28.2718 -3   -23.542 13.272 0.001794 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

Alpha Diversity

ds_patches %>%
  filter(disturbance == disturbance_input,
         patch_size == "S") %>%
  ggplot(
    aes(
      x = day,
      y = shannon,
      group = culture_ID,
      fill = culture_ID,
      color = culture_ID,
      linetype = eco_metaeco_type
    )
  ) +
  geom_line(stat = "summary", fun = "mean") +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Shannon Index",
    linetype = "",
    caption = "Vertical grey line: resource flow"
  ) +
  scale_linetype_discrete(labels = c(
    "Small isolated",
    "Small connected to large",
    "Small connected to small"
  ))  +
  scale_x_continuous(breaks = unique(ds_patches$day)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

ds_patches %>%
  filter(disturbance == disturbance_input,
         patch_size == "S") %>%
  ggplot(aes(
    x = day,
    y = shannon,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  labs(
    x = "Day",
    y = "Alpha Diversity (Shannon Index)",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_fill_manual(
    values = c(colour_isolated, colour_different_size, colour_same_size),
    labels = c(
      "Small isolated",
      "Small connected to large",
      "Small connected to small"
    )
  ) +
  scale_x_continuous(breaks = unique(ds_patches$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  )

p_effect_size = ds_patches_effect_size %>%
  filter(
    !time_point == 0,
    disturbance == disturbance_input,
    eco_metaeco_type %in% c("S (S_S)", "S (S_L)")
  ) %>%
  ggplot(aes(x = day,
             y = shannon_d,
             color = eco_metaeco_type)) +
  geom_point(position = position_dodge(0.5)) +
  geom_line(position = position_dodge(0.5)) +
  geom_errorbar(
    aes(ymin = shannon_d_lower,
        ymax = shannon_d_upper),
    width = .2,
    position = position_dodge(0.5)
  ) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted",
    color = "black",
    size = 0.7
  ) +
  labs(
    x = "Day",
    y = "Shannon Index Effect Size (Hedge's d)",
    color = ""
  ) +
  scale_color_discrete(labels = c("Small connected to large",
                                  "Small connnected to small")) +
  scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) 

p_effect_size

filtered_data = ds_patches_effect_size %>%
                         filter(time_point >= first_time_point_model,
                                time_point <= last_time_point_model,
                                disturbance == disturbance_input,
                                eco_metaeco_type %in% c("S (S_S)", "S (S_L)"))
full_model = lm(shannon_d ~                  
                  eco_metaeco_type + 
                  eco_metaeco_type : day,
                  data = filtered_data)

par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(shannon_d ~                  
                  1,
                  data = filtered_data)

anova = anova(full_model, null_model)
AIC = AIC(full_model, null_model)

anova
## Analysis of Variance Table
## 
## Model 1: shannon_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: shannon_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 11.232                           
## 2     11 22.014 -3   -10.782 2.5598  0.128
AIC
##            df      AIC
## full_model  5 43.26070
## null_model  2 45.33559
p_value = anova$`Pr(>F)`[2]

AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC =  AIC_full - AIC_null

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P = R2_full - R2_null

p_value = round(p_value, digits = 2)
R2_full = round(R2_full, digits = 2)
R2_P = round(R2_P, digits = 2)
deltaAIC = round(deltaAIC, digits = 2)

final_figure = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC, ", p = ", p_value, ", R2 (pach type) = ", R2_P))

final_figure

final_figure +
  theme_bw(base_size = presentation_size) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "none"
  ) +
  labs(caption = "")

ggsave(here("results", "presentations", "small_alpha.png"))

Species Densities

patch_size_input = "S"
Ble
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Ble,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Ble Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Ble_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Ble_d_lower,
            ymax = Ble_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Ble Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Ble_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Ble_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.54256 -0.20791 -0.03939  0.24344  0.67374 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  0.699145   0.480613   1.455    0.184
## eco_metaeco_typeS (S_S)     -0.527922   0.679689  -0.777    0.460
## eco_metaeco_typeS (S_L):day -0.037243   0.024963  -1.492    0.174
## eco_metaeco_typeS (S_S):day -0.008451   0.024963  -0.339    0.744
## 
## Residual standard error: 0.4177 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.2265, Adjusted R-squared:  -0.06363 
## F-statistic: 0.7807 on 3 and 8 DF,  p-value: 0.5372
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Ble_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 18.23815
## null_model  2 15.31937
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Ble_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Ble_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 1.3959                           
## 2     11 1.8046 -3  -0.40864 0.7807 0.5372
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Ble_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Ble_small

Cep
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Cep,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Cep Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Cep_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Cep_d_lower,
            ymax = Cep_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Cep Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Cep_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Cep_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69132 -0.31097 -0.06966  0.19031  1.02870 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -0.13153    0.65022  -0.202    0.845
## eco_metaeco_typeS (S_S)     -0.31176    0.91955  -0.339    0.743
## eco_metaeco_typeS (S_L):day  0.04444    0.03377   1.316    0.225
## eco_metaeco_typeS (S_S):day  0.02348    0.03377   0.695    0.507
## 
## Residual standard error: 0.5651 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.4548, Adjusted R-squared:  0.2503 
## F-statistic: 2.224 on 3 and 8 DF,  p-value: 0.1628
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Cep_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 25.49221
## null_model  2 26.77084
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Cep_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Cep_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 2.5550                           
## 2     11 4.6861 -3   -2.1311 2.2243 0.1628
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Cep_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Cep_small

Col
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Col,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Col Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Col_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Col_d_lower,
            ymax = Col_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Col Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Col_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Col_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8589 -0.4407  0.1131  0.3124  0.8175 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -0.17736    0.75190  -0.236    0.819
## eco_metaeco_typeS (S_S)     -0.01247    1.06335  -0.012    0.991
## eco_metaeco_typeS (S_L):day  0.05244    0.03905   1.343    0.216
## eco_metaeco_typeS (S_S):day  0.01035    0.03905   0.265    0.798
## 
## Residual standard error: 0.6535 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.4302, Adjusted R-squared:  0.2165 
## F-statistic: 2.013 on 3 and 8 DF,  p-value: 0.1908
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Col_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 28.97914
## null_model  2 29.72785
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Col_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Col_d ~ 1
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1      8 3.4165                          
## 2     11 5.9955 -3    -2.579 2.013 0.1908
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Col_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Col_small

Eug
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Eug,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Eug Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Eug_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Eug_d_lower,
            ymax = Eug_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Eug Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size
## Warning: Removed 7 rows containing missing values (geom_point).
## Warning: Removed 6 row(s) containing missing values (geom_path).

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Eug_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Eug_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##          7          9         13         14         15 
##  0.000e+00 -2.776e-17 -3.773e-01  7.545e-01 -3.773e-01 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -1.08953    2.06642  -0.527    0.691
## eco_metaeco_typeS (S_S)      0.64672    2.89790   0.223    0.860
## eco_metaeco_typeS (S_L):day  0.04559    0.16337   0.279    0.827
## eco_metaeco_typeS (S_S):day  0.00383    0.16337   0.023    0.985
## 
## Residual standard error: 0.9241 on 1 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.09765,    Adjusted R-squared:  -2.609 
## F-statistic: 0.03607 on 3 and 1 DF,  p-value: 0.9866
par(mfrow = c(2,3))
plot(full_model, which = 1:5)
## Warning: not plotting observations with leverage one:
##   1, 2
null_model = lm(Eug_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df       AIC
## full_model  5 15.353164
## null_model  2  9.866942
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Eug_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Eug_d ~ 1
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)
## 1      1 0.85402                           
## 2      4 0.94644 -3 -0.092422 0.0361 0.9866
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Eug_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Eug_small
## Warning: Removed 7 rows containing missing values (geom_point).
## Warning: Removed 6 row(s) containing missing values (geom_path).

Eup
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Eup,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Eup Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Eup_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Eup_d_lower,
            ymax = Eup_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Eup Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Eup_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Eup_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62632 -0.24352 -0.05568  0.37855  0.50876 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  0.187124   0.506564   0.369    0.721
## eco_metaeco_typeS (S_S)     -0.882553   0.716390  -1.232    0.253
## eco_metaeco_typeS (S_L):day  0.003398   0.026311   0.129    0.900
## eco_metaeco_typeS (S_S):day  0.023278   0.026311   0.885    0.402
## 
## Residual standard error: 0.4403 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.3875, Adjusted R-squared:  0.1578 
## F-statistic: 1.687 on 3 and 8 DF,  p-value: 0.2463
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Eup_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 19.50028
## null_model  2 19.38201
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Eup_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Eup_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 1.5507                           
## 2     11 2.5316 -3  -0.98091 1.6868 0.2463
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Eup_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Eup_small

Lox
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Lox,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Lox Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Lox_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Lox_d_lower,
            ymax = Lox_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Lox Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Lox_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Lox_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##        7        8        9       13       14       15       17       18 
## -0.27972  0.55944 -0.27972 -0.10389  0.63707 -0.78568  0.32822 -0.07571 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -0.79193    1.39775  -0.567    0.601
## eco_metaeco_typeS (S_S)      0.27675    1.57782   0.175    0.869
## eco_metaeco_typeS (S_L):day  0.06698    0.11239   0.596    0.583
## eco_metaeco_typeS (S_S):day  0.03388    0.03832   0.884    0.427
## 
## Residual standard error: 0.6358 on 4 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.2246, Adjusted R-squared:  -0.3569 
## F-statistic: 0.3863 on 3 and 4 DF,  p-value: 0.7697
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Lox_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 19.91137
## null_model  2 15.94665
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Lox_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Lox_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      4 1.6169                           
## 2      7 2.0853 -3  -0.46841 0.3863 0.7697
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Lox_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Lox_small
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).

Pau
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Pau,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Pau Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Pau_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Pau_d_lower,
            ymax = Pau_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Pau Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Pau_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Pau_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.51958 -0.08078  0.13608  0.30950  0.61788 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -0.758897   0.790054  -0.961    0.365
## eco_metaeco_typeS (S_S)     -0.004108   1.117306  -0.004    0.997
## eco_metaeco_typeS (S_L):day  0.062992   0.041036   1.535    0.163
## eco_metaeco_typeS (S_S):day  0.029295   0.041036   0.714    0.496
## 
## Residual standard error: 0.6867 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.3957, Adjusted R-squared:  0.1691 
## F-statistic: 1.746 on 3 and 8 DF,  p-value: 0.2349
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Pau_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 30.16710
## null_model  2 30.21141
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Pau_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Pau_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 3.7721                           
## 2     11 6.2421 -3     -2.47 1.7462 0.2349
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Pau_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Pau_small

Pca
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Pca,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Pca Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Pca_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Pca_d_lower,
            ymax = Pca_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Pca Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Pca_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Pca_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.18541 -0.18044  0.07589  0.27116  1.08334 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                  0.48406    0.84004   0.576   0.5803  
## eco_metaeco_typeS (S_S)      0.60926    1.18799   0.513   0.6219  
## eco_metaeco_typeS (S_L):day  0.08174    0.04363   1.873   0.0979 .
## eco_metaeco_typeS (S_S):day -0.05443    0.04363  -1.248   0.2475  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7301 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.7512, Adjusted R-squared:  0.6579 
## F-statistic: 8.053 on 3 and 8 DF,  p-value: 0.008432
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Pca_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 31.63938
## null_model  2 42.33425
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Pca_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Pca_d ~ 1
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
## 1      8  4.2644                                
## 2     11 17.1423 -3   -12.878 8.0529 0.008432 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Pca_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Pca_small

Spi
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Spi,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Spi Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Spi_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Spi_d_lower,
            ymax = Spi_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Spi Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Spi_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Spi_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38685 -0.15019 -0.01561  0.17460  0.53639 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                 -0.860297   0.367041  -2.344   0.0471 *
## eco_metaeco_typeS (S_S)      0.141441   0.519075   0.272   0.7921  
## eco_metaeco_typeS (S_L):day  0.038288   0.019064   2.008   0.0795 .
## eco_metaeco_typeS (S_S):day  0.005846   0.019064   0.307   0.7669  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.319 on 8 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.5531, Adjusted R-squared:  0.3855 
## F-statistic:   3.3 on 3 and 8 DF,  p-value: 0.07857
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Spi_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 11.76803
## null_model  2 15.43263
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Spi_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Spi_d ~ 1
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1      8 0.81413                              
## 2     11 1.82166 -3   -1.0075 3.3001 0.07857 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Spi_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Spi_small

Spi_te
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Spi_te,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Spi_te Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Spi_te_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Spi_te_d_lower,
            ymax = Spi_te_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Spi_te Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Spi_te_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Spi_te_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7785 -0.4592 -0.3043  0.5510  0.8552 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  0.09015    0.83646   0.108    0.917
## eco_metaeco_typeS (S_S)      0.17019    1.28491   0.132    0.898
## eco_metaeco_typeS (S_L):day  0.01087    0.04345   0.250    0.810
## eco_metaeco_typeS (S_S):day -0.02368    0.05747  -0.412    0.693
## 
## Residual standard error: 0.727 on 7 degrees of freedom
##   (7 observations deleted due to missingness)
## Multiple R-squared:  0.1333, Adjusted R-squared:  -0.2382 
## F-statistic: 0.3588 on 3 and 7 DF,  p-value: 0.7848
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(Spi_te_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 29.23044
## null_model  2 24.80390
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Spi_te_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Spi_te_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      7 3.6997                           
## 2     10 4.2686 -3  -0.56893 0.3588 0.7848
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Spi_te_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Spi_te_small
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

Tet
ds_patches %>%
      filter(
        disturbance == disturbance_input,
        patch_size == patch_size_input
      ) %>%
      ggplot(
        aes(
          x = day,
          y = Tet,
          group = interaction(day, eco_metaeco_type),
          fill = eco_metaeco_type
        )
      ) +
      geom_boxplot() +
      labs(
        x = "Day",
        y = "Tet Density (indiv/μl ?)",
        fill = "",
        caption = "Vertical grey line: resource flow"
      ) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      scale_fill_manual(
        values = c(colour_isolated, colour_different_size, colour_same_size),
        labels = c(
          "Small isolated",
          "Small connected to large",
          "Small connected to small"
        )
      )  +
      scale_x_continuous(breaks = unique(ds_patches$day)) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      )

p_effect_size = ds_patches_effect_size %>%
      filter(disturbance == disturbance_input,
             eco_metaeco_type %in% c("S (S_L)", "S (S_S)"),
             day > 0) %>%
      ggplot(aes(
        x = day,
        y = Tet_d,
        color = eco_metaeco_type
      )) +
      geom_point(position = position_dodge(0.5)) +
      geom_line(position = position_dodge(0.5)) +
      geom_errorbar(
        aes(ymin = Tet_d_lower,
            ymax = Tet_d_upper),
        width = .2,
        position = position_dodge(0.5)
      ) +
      labs(
        x = "Day",
        y = "Tet Density Effect Size (Hedge's d)",
        color = ""
      ) +
       scale_color_discrete(
         labels = c("Small connected to large",
                    "Small connnected to small")
       ) +
      scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
      theme_bw() +
      theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      ) +
      geom_vline(
        xintercept = resource_flow_days[1] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[2] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[3] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[4] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[5] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_vline(
        xintercept = resource_flow_days[6] + 0.7,
        linetype = "dotdash",
        color = "grey",
        size = 0.7
      ) +
      geom_hline(
        yintercept = 0,
        linetype = "dotted",
        color = "black",
        size = 0.7
      )

p_effect_size
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 7 row(s) containing missing values (geom_path).

filtered_data = ds_patches_effect_size %>%
  filter(patch_size == patch_size_input,
         disturbance == disturbance_input,
         time_point >= first_time_point_model,
         time_point <= last_time_point_model)
full_model = lm(Tet_d ~                  
                  eco_metaeco_type +
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = Tet_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##          7         10         11         13 
## -4.394e-02  1.758e-01 -1.318e-01  3.643e-17 
## 
## Coefficients: (1 not defined because of singularities)
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -1.14250    0.35427  -3.225    0.191
## eco_metaeco_typeS (S_S)      0.38370    0.41918   0.915    0.528
## eco_metaeco_typeS (S_L):day  0.07690    0.01903   4.041    0.154
## eco_metaeco_typeS (S_S):day       NA         NA      NA       NA
## 
## Residual standard error: 0.2241 on 1 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.9675, Adjusted R-squared:  0.9026 
## F-statistic:  14.9 on 2 and 1 DF,  p-value: 0.1802
par(mfrow = c(2,3))
plot(full_model, which = 1:5)
## Warning: not plotting observations with leverage one:
##   4
null_model = lm(Tet_d ~                  
                  1,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df       AIC
## full_model  4  1.839698
## null_model  2 11.548924
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)

anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: Tet_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: Tet_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      1 0.0502                           
## 2      3 1.5460 -2   -1.4958 14.897 0.1802
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

density_Tet_small = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

density_Tet_small
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 7 row(s) containing missing values (geom_path).

Large

Biomass

ds_patches %>%
  filter(disturbance == disturbance_input,
         patch_size == "L") %>%
  ggplot(
    aes(
      x = day,
      y = bioarea_per_volume,
      group = culture_ID,
      fill = culture_ID,
      color = culture_ID,
      linetype = eco_metaeco_type
    )
  ) +
  geom_line(stat = "summary", fun = "mean") +
  labs(
    x = "Day",
    y = "Patch Bioarea (µm²/μl)",
    linetype = "",
    caption = "Vertical grey line: resource flow"
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_linetype_discrete(labels = c(
    "Large isolated",
    "Large connected to large",
    "Large connected to small"
  ))  +
  scale_x_continuous(breaks = unique(ds_patches$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  )

ds_patches %>%
  filter(patch_size == "L",
         disturbance == disturbance_input) %>%
  ggplot(aes(
    x = day,
    y = bioarea_per_volume,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  labs(
    x = "Day",
    y = "Patch Bioarea (µm²/μl)",
    fill = ""
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_fill_manual(
    values = c(colour_isolated, colour_different_size, colour_same_size),
    labels = c(
      "Small isolated",
      "Small connected to large",
      "Small connected to small",
      caption = "Vertical grey line: resource flow"
    )
  ) +
  scale_x_continuous(breaks = unique(ds_patches$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) 

p_effect_size = ds_patches_effect_size %>%
  filter(
    !time_point == 0,
    disturbance == disturbance_input,
    eco_metaeco_type %in% c("L (L_L)", "L (S_L)")
  ) %>%
  ggplot(aes(x = day,
             y = bioarea_per_volume_d,
             color = eco_metaeco_type)) +
  geom_point(position = position_dodge(0.5)) +
  geom_line(position = position_dodge(0.5)) +
  geom_errorbar(
    aes(ymin = bioarea_per_volume_d_lower,
        ymax = bioarea_per_volume_d_upper),
    width = .2,
    position = position_dodge(0.5)
  ) +
  labs(
    x = "Day",
    y = "Bioarea Density Effect Size (Hedge's d)",
    color = ""
  ) +
  scale_color_discrete(labels = c("Small connected to large",
                                  "Small connnected to small")) +
  scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted",
    color = "black",
    size = 0.7
  )

p_effect_size

filtered_data = ds_patches_effect_size %>%
                         filter(time_point >= first_time_point_model,
                                time_point <= last_time_point_model,
                                disturbance == disturbance_input,
                                eco_metaeco_type %in% c("L (L_L)", "L (S_L)"))
full_model = lm(bioarea_per_volume_d ~                  
                  eco_metaeco_type + 
                  day : eco_metaeco_type,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = bioarea_per_volume_d ~ eco_metaeco_type + day:eco_metaeco_type, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70065 -0.36021  0.05607  0.28222  0.87208 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                  0.76759    0.67687   1.134   0.2896  
## eco_metaeco_typeL (S_L)     -0.17902    0.95723  -0.187   0.8563  
## eco_metaeco_typeL (L_L):day -0.07419    0.03516  -2.110   0.0678 .
## eco_metaeco_typeL (S_L):day -0.07262    0.03516  -2.066   0.0727 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5883 on 8 degrees of freedom
## Multiple R-squared:  0.5271, Adjusted R-squared:  0.3497 
## F-statistic: 2.972 on 3 and 8 DF,  p-value: 0.09688
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(bioarea_per_volume_d ~                  
                  day,
                  data = filtered_data)

AIC = AIC(full_model, null_model)
AIC
##            df      AIC
## full_model  5 26.45599
## null_model  3 22.74915
AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC = AIC_full - AIC_null
deltaAIC = round(deltaAIC, digits = 2)


anova = anova(full_model, null_model)
anova
## Analysis of Variance Table
## 
## Model 1: bioarea_per_volume_d ~ eco_metaeco_type + day:eco_metaeco_type
## Model 2: bioarea_per_volume_d ~ day
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 2.7687                           
## 2     10 2.8371 -2 -0.068471 0.0989 0.9069
p_value = anova$`Pr(>F)`[2]
p_value = round(p_value, digits = 3)

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P =  R2_full - R2_null
R2_P = round(R2_P, digits = 2)

large_biomass = p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC , ", p = ", p_value, ", Adjusted R2 = ", R2_P))

large_biomass

Alpha Diversity

ds_patches %>%
  filter(disturbance == disturbance_input,
         patch_size == "L") %>%
  ggplot(
    aes(
      x = day,
      y = shannon,
      group = culture_ID,
      fill = culture_ID,
      color = culture_ID,
      linetype = eco_metaeco_type
    )
  ) +
  geom_line(stat = "summary", fun = "mean") +
  labs(
    x = "Day",
    y = "Shannon Index",
    linetype = "",
    caption = "Vertical grey line: resource flow"
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_linetype_discrete(labels = c(
    "Large isolated",
    "Large connected to large",
    "Large connected to small"
  ))  +
  scale_x_continuous(breaks = unique(ds_patches$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  )

ds_patches %>%
  filter(disturbance == disturbance_input,
         patch_size == "L") %>%
  ggplot(aes(
    x = day,
    y = shannon,
    group = interaction(day, eco_metaeco_type),
    fill = eco_metaeco_type
  )) +
  geom_boxplot() +
  labs(
    x = "Day",
    y = "Shannon Index",
    fill = "",
    caption = "Vertical grey line: resource flow"
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  ) +
  scale_fill_manual(
    values = c(colour_isolated, colour_same_size, colour_different_size),
    labels = c(
      "Small isolated",
      "Large connected to large",
      "Large connected to small"
    )
  ) +
  scale_x_continuous(breaks = unique(ds_patches$day)) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  )

p_effect_size = ds_patches_effect_size %>%
  filter(
    !time_point == 0,
    disturbance == disturbance_input,
    eco_metaeco_type == "L (L_L)" |
      eco_metaeco_type == "L (S_L)"
  ) %>%
  ggplot(aes(x = day,
             y = shannon_d,
             color = eco_metaeco_type)) +
  geom_point(position = position_dodge(0.5)) +
  geom_line(position = position_dodge(0.5)) +
  geom_errorbar(
    aes(ymin = shannon_d_lower,
        ymax = shannon_d_upper),
    width = .2,
    position = position_dodge(0.5)
  ) +
  geom_vline(
    xintercept = resource_flow_days[1] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[2] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[3] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[4] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[5] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_vline(
    xintercept = resource_flow_days[6] + 0.7,
    linetype = "dotdash",
    color = "grey",
    size = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted",
    color = "black",
    size = 0.7
  ) +
  scale_x_continuous(breaks = unique(ds_patches_effect_size$day)) +
  labs(
    x = "Day",
    y = "Shannon Index Effect Size (Hedge's d)",
    color = ""
  ) +
  scale_color_discrete(labels = c("Large connected to large",
                                  "Large connected to small")) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

p_effect_size

filtered_data = ds_patches_effect_size %>%
                         filter(time_point >= first_time_point_model,
                                time_point <= last_time_point_model,
                                disturbance == disturbance_input,
                                eco_metaeco_type %in% c("L (L_L)", "L (S_L)"))
full_model = lm(shannon_d ~                  
                  eco_metaeco_type + 
                  eco_metaeco_type : day,
                  data = filtered_data)

summary(full_model)
## 
## Call:
## lm(formula = shannon_d ~ eco_metaeco_type + eco_metaeco_type:day, 
##     data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93436 -0.41622 -0.07603  0.52024  0.92419 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 -1.21097    0.90758  -1.334    0.219
## eco_metaeco_typeL (S_L)     -0.28707    1.28351  -0.224    0.829
## eco_metaeco_typeL (L_L):day  0.04688    0.04714   0.995    0.349
## eco_metaeco_typeL (S_L):day  0.05848    0.04714   1.240    0.250
## 
## Residual standard error: 0.7888 on 8 degrees of freedom
## Multiple R-squared:  0.2422, Adjusted R-squared:  -0.04192 
## F-statistic: 0.8525 on 3 and 8 DF,  p-value: 0.5034
par(mfrow = c(2,3))
plot(full_model, which = 1:5)

null_model = lm(shannon_d ~                  
                  1,
                  data = filtered_data)

anova = anova(full_model, null_model)
AIC = AIC(full_model, null_model)

anova
## Analysis of Variance Table
## 
## Model 1: shannon_d ~ eco_metaeco_type + eco_metaeco_type:day
## Model 2: shannon_d ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 4.9778                           
## 2     11 6.5691 -3   -1.5913 0.8525 0.5034
AIC
##            df      AIC
## full_model  5 33.49542
## null_model  2 30.82413
p_value = anova$`Pr(>F)`[2]

AIC_full = AIC$AIC[1]
AIC_null = AIC$AIC[2]
deltaAIC =  AIC_full - AIC_null

R2_full = glance(full_model)$r.squared
R2_null = glance(null_model)$r.squared
R2_P = R2_full - R2_null

p_value = round(p_value, digits = 2)
R2_full = round(R2_full, digits = 2)
R2_P = round(R2_P, digits = 2)
deltaAIC = round(deltaAIC, digits = 2)

p_effect_size +
  labs(title = paste0("ΔAIC = ", deltaAIC, ", p = ", p_value, ", R2 (pach type) = ", R2_P))

Evaporation

Effects of evaporation

We want to know if there was a systematic bias in the evaporation of different treatments (disturbance, patch size) and whether evaporation changed across time. My expectation would be that we would see a difference among the exchanges 2,3 and the exchanges 4,5,6. This is because in exchange 2,3 cultures were microwaved in 15 tubes for 3 minutes and in exchange 4,5,6 cultures were microwaved in 4 tubes for only 1 minute.

Tidy

#Columns: exchange & evaporation
ds_for_evaporation = gather(ds_for_evaporation, 
                            key = exchange, 
                            value = evaporation, 
                            water_add_after_t2:water_add_after_t6)
ds_for_evaporation[ds_for_evaporation == "water_add_after_t2"] = "2"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t3"] = "3"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t4"] = "4"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t5"] = "5"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t6"] = "6"
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] / 2 #This is because exchange contained the topping up of two exchanges
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] + 2 #We need to add 2 ml to the evaporation that happened at the exchange events 1 and 2. This is because we already added 1 ml of water at exchange 1 and 1 ml of water at exchange 2. 

#Column: nr_of_tubes_in_rack
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 1] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 2] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 3] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 4] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 5] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 6] = 4

Plot

ds_for_evaporation %>%
  filter(disturbance == disturbance) %>%
  ggplot(aes(x = as.character(nr_of_tubes_in_rack),
             y = evaporation)) + 
  geom_boxplot() + 
  labs(x = "Number of tubes in rack", 
       y = "Evaporation (ml)")

ds_for_evaporation %>%
  filter(disturbance == disturbance) %>%
  ggplot(aes(x = as.character(patch_size),
             y = evaporation)) + 
  geom_boxplot() + 
  labs(x = "Patch size", 
       y = "Evaporation (ml)")

ds_for_evaporation %>%
  filter(disturbance == disturbance) %>%
  ggplot(aes(x = as.character(day),
             y = evaporation)) + 
  geom_boxplot() + 
  labs(x = "Day", 
       y = "Evaporation (ml)")

ds_for_evaporation %>%
  filter(disturbance == disturbance) %>%
  ggplot(aes(x = disturbance,
             y = evaporation)) + 
  geom_boxplot() + 
  labs(x = "Disturbance", 
       y = "Evaporation (ml)")

It seems like there is no real difference across time, disturbance, or patch type. However, we could also run a mixed effect model to show that they do not.

Mixed effect model

It gives me the following error:

  • Error in fn(nM$xeval()) : Downdated VtV is not positive definite
mixed.model = lmer(evaporation  ~ 
                     patch_size * disturbance  * exchange + 
                     (exchange | culture_ID), 
                   data = ds_for_evaporation,
                   REML = FALSE, 
                   control = lmerControl (optimizer = "Nelder_Mead"))

null.model = lm(evaporation ~
                  1, 
                data = ds_for_evaporation)

anova(mixed.model, null.model)

Evaporation tests

Evaporation when microwaving 15 falcon tubes at the time

evaporation.test = read.csv(here("data", "evaporation_test","evaporation_test_right.csv"), header = TRUE)

evaporation.test %>%
  ggplot(aes (x = as.character(water_pipetted),
                y = weight_water_evaporated,
                group = interaction(water_pipetted, as.character(rack)),
                fill = as.character(rack))) +
  geom_boxplot() +
  labs(x = "Water volume (ml)" , 
       y = "Evaporation (g)", 
       fill = "Rack replicate")

Evaporation when microwaving 5 tubes with 10 filled or empty tubes

evaporation.test = read.csv(here("data", "evaporation_test", "evaporation_test_fill_nofill.csv"), header = TRUE)

evaporation.test %>%
  ggplot(aes (x = all_tubes_water,
              y = weight_water_evaporated)) +
  geom_boxplot() +
  labs(x = "Water in the other 10 tubes" , 
  y = "Evaporation (g)", 
  caption = "When all tubes were filled, they were filled with 6.75 ml of deionised water.")

Video analysis script

rm(list = ls())
setwd("/media/mendel-himself/ID_061_Ema2/PatchSizePilot/training")

# load package
# library(devtools)
# install_github("femoerman/bemovi", ref="master")
library(bemovi)
library(parallel)
library(doParallel)
library(foreach)

#Define memory to be allocated
memory.alloc <- 240000 #-needs_to_be_specified
memory.per.identifier <- 40000 #-needs_to_be_specified
memory.per.linker <- 5000 #-needs_to_be_specified
memory.per.overlay <- 60000 #-needs_to_be_specified

# UNIX
# set paths to tools folder and particle linker
tools.path <-
  "/home/mendel-himself/bemovi_tools/" #-needs_to_be_specified
to.particlelinker <- tools.path

# directories and file names
to.data <- paste(getwd(), "/", sep = "")
video.description.folder <- "0_video_description/"
video.description.file <- "video_description.txt"
raw.video.folder <- "1_raw/"
raw.avi.folder <- "1a_raw_avi/"
metadata.folder <- "1b_raw_meta/"
particle.data.folder <- "2_particle_data/"
trajectory.data.folder <- "3_trajectory_data/"
temp.overlay.folder <- "4a_temp_overlays/"
overlay.folder <- "4_overlays/"
merged.data.folder <- "5_merged_data/"
ijmacs.folder <- "ijmacs/"


######################################################################
# VIDEO PARAMETERS

# video frame rate (in frames per second)
fps <- 25 #-needs_to_be_specified

# length of video (in frames)
total_frames <- 125 #-needs_to_be_specified

#Dimensions of the videos in pixels
width = 2048 #-needs_to_be_specified
height = 2048 #-needs_to_be_specified

# measured volume (in microliter) #-needs_to_be_specified
measured_volume <-
  34.4 # for Leica M205 C with 1.6 fold magnification, sample height 0.5 mm and Hamamatsu Orca Flash 4
#measured_volume <- 14.9 # for Nikon SMZ1500 with 2 fold magnification, sample height 0.5 mm and Canon 5D Mark III

# size of a pixel (in micrometer) #-needs_to_be_specified
pixel_to_scale <-
  4.05 # for Leica M205 C with 1.6 fold magnification, sample height 0.5 mm and Hamamatsu Orca Flash 4
#pixel_to_scale <- 3.79 # for Nikon SMZ1500 with 2 fold magnification, sample height 0.5 mm and Canon 5D Mark III

# specify video file format (one of "avi","cxd","mov","tiff")
# bemovi only works with avi and cxd. other formats are reformated to avi below
video.format <- "cxd" #-needs_to_be_specified

# setup
difference.lag <- 10
thresholds <- c(13, 255) # don't change the second value
# thresholds <- c(50,255)

# MORE PARAMETERS (USUALLY NOT CHANGED)
######################################################################
# FILTERING PARAMETERS
# optimized for Perfex Pro 10 stereomicrocope with Perfex SC38800 (IDS UI-3880LE-M-GL) camera
# tested stereomicroscopes: Perfex Pro 10, Nikon SMZ1500, Leica M205 C
# tested cameras: Perfex SC38800, Canon 5D Mark III, Hamamatsu Orca Flash 4
# tested species: Tet, Col, Pau, Pca, Eug, Chi, Ble, Ceph, Lox, Spi

# min and max size: area in pixels
particle_min_size <- 10
particle_max_size <- 1000

# number of adjacent frames to be considered for linking particles
trajectory_link_range <- 3
# maximum distance a particle can move between two frames
trajectory_displacement <- 16

# these values are in the units defined by the parameters above: fps (seconds), measured_volume (microliters) and pixel_to_scale (micometers)
filter_min_net_disp <- 25
filter_min_duration <- 1
filter_detection_freq <- 0.1
filter_median_step_length <- 3

######################################################################
# VIDEO ANALYSIS

#Check if all tools are installed, and if not install them
check_tools_folder(tools.path)

#Ensure computer has permission to run bftools
system(paste0("chmod a+x ", tools.path, "bftools/bf.sh"))
system(paste0("chmod a+x ", tools.path, "bftools/bfconvert"))
system(paste0("chmod a+x ", tools.path, "bftools/showinf"))

# Convert files to compressed avi (takes approx. 2.25 minutes per video)
convert_to_avi(
  to.data,
  raw.video.folder,
  raw.avi.folder,
  metadata.folder,
  tools.path,
  fps,
  video.format
)


# TESTING

# check file format and naming
# check_video_file_names(to.data,raw.avi.folder,video.description.folder,video.description.file)

# check whether the thresholds make sense (set "dark backgroud" and "red")
#check_threshold_values(to.data, raw.avi.folder, ijmacs.folder, 2, difference.lag, thresholds, tools.path,  memory.alloc)

# identify particles
locate_and_measure_particles(
  to.data,
  raw.avi.folder,
  particle.data.folder,
  difference.lag,
  min_size = particle_min_size,
  max_size = particle_max_size,
  thresholds = thresholds,
  tools.path,
  memory = memory.alloc,
  memory.per.identifier = memory.per.identifier,
  max.cores = detectCores() - 1
)

# link the particles
link_particles(
  to.data,
  particle.data.folder,
  trajectory.data.folder,
  linkrange = trajectory_link_range,
  disp = trajectory_displacement,
  start_vid = 1,
  memory = memory.alloc,
  memory_per_linkerProcess = memory.per.linker,
  raw.avi.folder,
  max.cores = detectCores() - 1,
  max_time = 1
)

# merge info from description file and data
merge_data(
  to.data,
  particle.data.folder,
  trajectory.data.folder,
  video.description.folder,
  video.description.file,
  merged.data.folder
)

# load the merged data
load(paste0(to.data, merged.data.folder, "Master.RData"))

# filter data: minimum net displacement, their duration, the detection frequency and the median step length
trajectory.data.filtered <-
  filter_data(
    trajectory.data,
    filter_min_net_disp,
    filter_min_duration,
    filter_detection_freq,
    filter_median_step_length
  )

# summarize trajectory data to individual-based data
morph_mvt <-
  summarize_trajectories(
    trajectory.data.filtered,
    calculate.median = F,
    write = T,
    to.data,
    merged.data.folder
  )

# get sample level info
summarize_populations(
  trajectory.data.filtered,
  morph_mvt,
  write = T,
  to.data,
  merged.data.folder,
  video.description.folder,
  video.description.file,
  total_frames
)

# create overlays for validation
create.subtitle.overlays(
  to.data,
  traj.data = trajectory.data.filtered,
  raw.video.folder,
  raw.avi.folder,
  temp.overlay.folder,
  overlay.folder,
  fps,
  vid.length = total_frames / fps,
  width,
  height,
  tools.path = tools.path,
  overlay.type = "number",
  video.format
)

# Create overlays (old method)
create_overlays(
  traj.data = trajectory.data.filtered,
  to.data = to.data,
  merged.data.folder = merged.data.folder,
  raw.video.folder = raw.avi.folder,
  temp.overlay.folder = "4a_temp_overlays_old/",
  overlay.folder = "4_overlays_old/",
  width = width,
  height = height,
  difference.lag = difference.lag,
  type = "traj",
  predict_spec = F,
  contrast.enhancement = 1,
  IJ.path = "/home/mendel-himself/bemovi_tools",
  memory = memory.alloc,
  max.cores = detectCores() - 1,
  memory.per.overlay = memory.per.overlay
)

########################################################################
# some cleaning up
#system("rm -r 2_particle_data")
#system("rm -r 3_trajectory_data")
#system("rm -r 4a_temp_overlays")
system("rm -r ijmacs")
########################################################################

Here I want to check whether I can distinguish between Tet and Chi (food of Ble). To do so, let’s look at the traits of Ble and Chi.

load(here("data", "morphology", "training.RData"))
training = morph_mvt
rm(morph_mvt)
training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_grey)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_area)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_perimeter)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_major)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_minor)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_ar)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = mean_turning)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = duration)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = N_frames)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = max_net)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = net_disp)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = net_speed)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = gross_disp)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = gross_speed)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = max_step)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = min_step)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = max_gross_speed)) +
  geom_boxplot()

training %>%
  filter(species == "Ble" | species == "Tet") %>%
  ggplot(aes(x = species,
             y = min_gross_speed)) +
  geom_boxplot()

It seems like there’s no great way of telling them apart.

Figures

## Warning: Removed 7 rows containing missing values (geom_point).
## Warning: Removed 6 row(s) containing missing values (geom_path).

## Warning: Removed 7 rows containing missing values (geom_point).
## Removed 6 row(s) containing missing values (geom_path).

## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).

## Warning: Removed 4 rows containing missing values (geom_point).
## Removed 3 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).
## Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 7 row(s) containing missing values (geom_path).

## Warning: Removed 9 rows containing missing values (geom_point).
## Removed 7 row(s) containing missing values (geom_path).

Running time

## Time difference of 1.357676 mins

Bibliography